/** 
 * @main Plik z definicjami klasy GenLosowy.
 * @file    losowy.cpp
 * @author  Grzegorz Sejnota
 *
 * @created on 5 styczen 2009, 19:53
 * @version 0.0001
 */

#include "losowy.h"
#include <iostream>
#include <cmath>
#include <ctime>
using namespace std;

template <typename T>
void wypiszVector (vector<T> wektor){
    for(int i=0; i < wektor.size(); i++){
        cout<< endl <<"\t" << wektor[i];
    }
}

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//++++++++++++++ POCZĄTEK [klasy GenLosowy] ++++++++++++++++++++++++++++++
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//======= konstruktor ==========================================================
/**
 * Konstruktor prywatny, wykorzystywany w metodzie pobierzInstancje().
 *
 *      @param dlugoscCiagu - dlugosc ciagu liczb losowych o wariancji 1
 */
GenLosowy::GenLosowy(int dlugoscCiagu) {
    double liczbaLosowa;
    numerProbki = 0;
    srand(time(NULL));
    for (int i = 0; i < dlugoscCiagu; i++) {
        liczbaLosowa = static_cast<double>( rand() ) / (static_cast<double> (RAND_MAX) + 1.0) 
                *2.0 - 1.0;

        ciagLosowy.push_back(liczbaLosowa);
    }
    //cout<<"ciag:"<<endl;
   //wypiszVector(ciagLosowy);
    //cout<<"wariancja ciagu:"<<endl;
    double var = obliczWariancje(ciagLosowy);
}

//======= destruktor ===========================================================
/**
 * Destruktor
 */
GenLosowy :: ~GenLosowy(){
}

//======= konstruktor kopiujacy ================================================
/**
 * Pusty konstruktor kopiujacy - stworzony by nie zostal utworzony konstruktor
 * kopiujacy przez kompilator.
 */
GenLosowy :: GenLosowy(const GenLosowy &){
}

//======= operator przypisania =================================================
/**
 * Pusty operator przypisania.
 */
/*
GenLosowy & GenLosowy :: operator=(const GenLosowy&){
}*/

//======= pobierzInstancje =====================================================
/**
 * Metoda tworzaca obiekt klasy GenLosowy lub zwracajaca referencje do
 * jedynego istniejacego jej obiektu - singletonu.
 *      @return referencje do jedynego obiektu klasy PlikKonfiguracyjny
 */
GenLosowy& GenLosowy :: pobierzInstancje(){
          //Inicjalizacja statycznego obiektu.
          //Obiekt zostanie utworzony przy pierwszym wywołaniu tej metody
          //i tylko wtedy nastąpi inicjalizacja przy pomocy konstruktora.
          //Każde następne wywołanie jedynie zwróci referencję tego obiektu.
    static GenLosowy instance = GenLosowy(10000);
    return instance;
}

//======= obliczWariancje ======================================================
/**
 * Metoda obliczajaca wartosc wariancji ciagu liczbowego.
 *      @param ciag - vector ciagu liczb
 *      @return wariancja ciagu liczb
 */
double GenLosowy :: obliczWariancje(const vector<double> & ciag){
    // ---------- obliczanie sredniej ------------------------------------------
    double srednia = 0.0;
    for(int i=0; i<ciag.size(); i++){
        srednia += ciag.at(i);
    }
    srednia /= ciag.size();
    //cout<<"srednia = "<<srednia<<endl;
    
    // ---------- obliczanie wariancji -----------------------------------------
    double wariancja = 0.0;
    for(int i=0; i<ciag.size(); i++){
        wariancja += pow((ciag.at(i)-srednia), 2);
    }
    wariancja /= ciag.size();
    //cout<<"dlugosc = "<<ciag.size()<<"wariancja = "<<wariancja<<endl;
    return wariancja;
}

//======= filtrujArma ==========================================================
/**
 * Metoda filtrujaca ciagLosowy przez filtr ARMA. Filtr ARMA (auto regresive 
 * moving average) jest zdefiniowany następująco:<br>
 *      y = C(z^-1)/A(z^-1)* e
 *
 *      @param wielomianA - wielomian w mianowniku transmitancji ARMA
 *      @param wielomianC - wielomian w liczniku transmitancji ARMA
 *      @return wariancja ciagu liczb
 */
vector<double> GenLosowy :: filtrujArma(Wielomian wielomianA, 
                                              Wielomian wielomianC){
    vector<double> ciagWyjsciowy(ciagLosowy.size(), 0.0);
    double y = 0.0;
    //cout<<endl<<"filtruje ARMA"<<endl;
    for(int j=0; j<ciagLosowy.size(); j++){
        double y = 0.0;
        //cout<<"j = "<<j<<endl;
        //-------- wplyw wyjsc -------------------------------------------------
        for(int i=0; i<wielomianA.pobierzStopien(); i++){
            //cout<<"indeks = "<<j-wielomianA.pobierzStopien()+i<<"  i = "<<i<<endl;
            if( (j-wielomianA.pobierzStopien()+i) >= 0 ) //aby nie odwolac sie do ujemnego indeksu
                y -= wielomianA[i]*ciagWyjsciowy.at(j-wielomianA.pobierzStopien()+i);
            //cout<<"ji="<<j<<i<<"   y="<<y<<endl;
        }
        //-------- wplyw wejsc -------------------------------------------------
        for(int i=0; i<=wielomianC.pobierzStopien(); i++){
            //cout<<"indeks = "<<j-wielomianC.pobierzStopien()+i<<endl;
            if( (j-wielomianC.pobierzStopien()+i) >= 0 ) //aby nie odwolac sie do ujemnego indeksu
                y += wielomianC[i]*ciagLosowy.at(j-wielomianC.pobierzStopien()+i);
            //cout<<"ji="<<j<<i<<"   y="<<y<<endl;
        }
        ciagWyjsciowy.at(j)= y ; //dodanie obliczonego wyjscia do wektora wyjsciowego
    }
    //wypiszVector(ciagWyjsciowy);
    
    
    /*for(int j=0; j<ciagWyjsciowy.size(); j++){
        ciagWariancja1.at(j) = wsp*ciagWariancja1.at(j);
    }
    obliczWariancje(ciagWyjsciowy);*/
    return ciagWyjsciowy;
}

//======= filtrujArma ==========================================================
/**
 * Metoda filtrujaca podany ciagWejsciowy przez filtr ARMA. Filtr ARMA (auto 
 * regresive moving average) jest zdefiniowany następująco:<br>
 *      y = C(z^-1)/A(z^-1)* e
 *
 *      @param wielomianA - wielomian w mianowniku transmitancji ARMA
 *      @param wielomianC - wielomian w liczniku transmitancji ARMA
 *      @param ciagWejsciowy - ciag 
 *      @return wariancja ciagu liczb
 */
vector<double> GenLosowy :: filtrujArma(Wielomian wielomianA, 
                                              Wielomian wielomianC,
                                              const vector<double> & ciagWejsciowy){
    vector<double> ciagWyjsciowy(ciagWejsciowy.size(), 0.0);
    
    //cout<<endl<<"filtruje ARMA"<<endl;
    for(int j=0; j<ciagWejsciowy.size(); j++){
        double y = 0.0;
        //cout<<"j = "<<j<<endl;
        //-------- wplyw wyjsc -------------------------------------------------
        for(int i=0; i<wielomianA.pobierzStopien(); i++){
            //cout<<"indeks = "<<j-wielomianA.pobierzStopien()+i<<"  i = "<<i<<endl;
            if( (j-wielomianA.pobierzStopien()+i) >= 0 ) //aby nie odwolac sie do ujemnego indeksu
                y -= wielomianA[i]*ciagWyjsciowy.at(j-wielomianA.pobierzStopien()+i);
            //cout<<"ji="<<j<<i<<"   y="<<y<<endl;
        }
        //-------- wplyw wejsc -------------------------------------------------
        for(int i=0; i<=wielomianC.pobierzStopien(); i++){
            //cout<<"indeks = "<<j-wielomianC.pobierzStopien()+i<<endl;
            if( (j-wielomianC.pobierzStopien()+i) >= 0 ) //aby nie odwolac sie do ujemnego indeksu
                y += wielomianC[i]*ciagWejsciowy.at(j-wielomianC.pobierzStopien()+i);
            //cout<<"ji="<<j<<i<<"   y="<<y<<endl;
        }
        ciagWyjsciowy.at(j)= y ; //dodanie obliczonego wyjscia do wektora wyjsciowego
    }
    //wypiszVector(ciagWyjsciowy);
    //double varPrzefiltrowanego = obliczWariancje(ciagWyjsciowy);
    
    /*for(int j=0; j<ciagWyjsciowy.size(); j++){
        ciagWariancja1.at(j) = wsp*ciagWariancja1.at(j);
    }
    obliczWariancje(ciagWyjsciowy);*/
    return ciagWyjsciowy;
}

//======= obliczWspPrzeskalowaniaCiagu =========================================
/**
 * Metoda obliczajaca wspolczynnik przeskalowania ciagLosowy, ktory pozwala
 * uzyskac wyjscie z toru zaklocenia ObiektARMAX o zadanej wariancji.
 *
 *      @param wariancjaZadana - zadana wariancja wyjscia toru zaklocenia
 *      @param wielomianA - wielomian w mianowniku transmitancji ARMA
 *      @param wielomianC - wielomian w liczniku transmitancji ARMA
 *      @return wariancja ciagu liczb
 */
double GenLosowy :: obliczWspPrzeskalowaniaCiagu(double wariancjaZadana, 
                                                       Wielomian wielomianA, 
                                                       Wielomian wielomianC){
    if (wariancjaZadana != 0.0){
        //-------- filtracja ciaguLosowego -----------------------------------------
        vector<double> ciagPrzefiltrowany = filtrujArma(wielomianA, wielomianC);
         //-------- obliczenie wspPrzeskalowania -----------------------------------
        double wariancjaPrzefiltrowanego = obliczWariancje(ciagPrzefiltrowany);
        double wspPrzeskalowania = sqrt(wariancjaZadana / wariancjaPrzefiltrowanego);

 /*
        cout<<endl<<"WSP PRZ = "<<wspPrzeskalowania<<endl;
        vector<double> ciagPrzeskalowany(ciagLosowy.size(), 0.0);
        for(int j=0; j<ciagLosowy.size(); j++){
            ciagPrzeskalowany.at(j) = wspPrzeskalowania*ciagLosowy.at(j);
        }
        cout<<endl<<"wariancja przeskalowanego i przefiltrowanego!"<<endl;
        obliczWariancje( filtrujArma(wielomianA, wielomianC, ciagPrzeskalowany) );
        //wypiszVector(filtrujArma(wielomianA, wielomianC, ciagPrzeskalowany));
 */       
        return wspPrzeskalowania;
    }else{
        return 0.0;
    }
}

//======= pobierzLiczbeLosowa =========================================
/**
 * Metoda zwracajaca wartosc probki z ciagu ciagLosowy na numerze numerProbki.
 */
double GenLosowy :: pobierzLiczbeLosowa(){
    double losowa = ciagLosowy.at(numerProbki);
    numerProbki++;
    if (numerProbki >= ciagLosowy.size() )
        numerProbki = 0;
    return losowa;
}

//------------------------------------------------------------------------------
//-------------- KONIEC [klasy GenLosowy] --------------------------------
//------------------------------------------------------------------------------

